Paso 5 (1)- Trayectorias de hospitalización y mortalidad con foco en condiciones vinculadas a trastornos de salud mental y consumo de sustancias posterior a un primer ingreso por alguno de estos trastornos, en usuarios/as jóvenes y adultos emergentes de población general y pertenecientes a pueblos originarios, 2018-2021, Chile
Representar las mejores opciones de agrupamiento, junto con su relación con otras variables. En esta etapa, analizamos la asociación con covariables de la solución con mejores índices de calidad para la resolución trimestral.
Autor/a
Andrés González Santa Cruz
Fecha de publicación
13 de may, 2025
Configurar
Código
# remover objetos y memoria utilizadarm(list=ls());gc()
used (Mb) gc trigger (Mb) max used (Mb)
Ncells 600033 32.1 1288752 68.9 756906 40.5
Vcells 1175765 9.0 8388608 64.0 1876026 14.4
#elegir repositorioif(Sys.info()["sysname"]=="Windows"){options(repos =c(CRAN ="https://cran.dcc.uchile.cl/"))}options(install.packages.check.source ="yes") # Chequea la fuente de los paquetes#borrar caché#system("fc-cache -f -v")if(!require(pacman)){install.packages("pacman");require(pacman)}pacman::p_unlock(lib.loc =.libPaths()) #para no tener problemas reinstalando paquetesif(Sys.info()["sysname"]=="Windows"){if (getRversion() !="4.4.1") { stop("Requiere versión de R 4.4.1. Actual: ", getRversion()) }}cat("quarto version: "); system("quarto --version")
quarto version:
[1] 0
Código
if(!require(job)){install.packages("job");require(job)}if(!require(kableExtra)){install.packages("kableExtra");require(kableExtra)}if(!require(tidyverse)){install.packages("tidyverse");require(tidyverse)}if(!require(cluster)){install.packages("cluster"); require(cluster)}if(!require(WeightedCluster)){install.packages("WeightedCluster"); require(WeightedCluster)}if(!require(devtools)){install.packages("devtools"); require(devtools)}if(!require(TraMineR)){install.packages("TraMineR"); require(TraMineR)}if(!require(TraMineRextras)){install.packages("TraMineRextras"); require(TraMineRextras)}if(!require(NbClust)){install.packages("NbClust"); require(NbClust)}if(!require(haven)){install.packages("haven"); require(haven)}if(!require(ggseqplot)){install.packages("ggseqplot"); require(ggseqplot)}if(!require(grid)){install.packages("grid"); require(grid)}if(!require(gridExtra)){install.packages("gridExtra"); require(gridExtra)}if(!require(Tmisc)){install.packages("Tmisc"); require(Tmisc)}if(!require(factoextra)){install.packages("factoextra"); require(factoextra)}if(!require(stargazer)){install.packages("stargazer"); require(stargazer)}if(!require(gtsummary)){install.packages("gtsummary"); require(gtsummary)}if(!require(lmtest)){install.packages("lmtest"); require(lmtest)}if(!require(emmeans)){install.packages("emmeans"); require(emmeans)}if(!require(fpp2)){install.packages("fpp2"); require(fpp2)}if(!require(purrr)){install.packages("purrr"); require(purrr)}if(!require(forecast)){install.packages("forecast"); require(forecast)}if(!require(magrittr)){install.packages("magrittr"); require(magrittr)}if(!require(foreach)){install.packages("foreach"); require(foreach)}if(!require(doParallel)){install.packages("doParallel"); require(doParallel)}if(!require(progressr)){install.packages("progressr"); require(progressr)}if(!require(chisq.posthoc.test)){devtools::install_github("ebbertd/chisq.posthoc.test")}if(!require(rstatix)){install.packages("rstatix"); require(rstatix)}if(!require(rio)){install.packages("rio"); require(rio)}if(!require(cowplot)){install.packages("cowplot"); require(cowplot)}if(!require(DiagrammeR)){install.packages("DiagrammeR"); require(DiagrammeR)}if(!require(DiagrammeRsvg)){install.packages("DiagrammeRsvg"); require(DiagrammeRsvg)}if(!require(rsvg)){install.packages("rsvg"); require(rsvg)}if(!require(survminer)){install.packages("survminer"); require(survminer)}if(!require(epitools)){install.packages("epitools"); require(epitools)}seq_mean_t_dos_grupos <-function(bd =NULL, group1, group2) {# Agrupar por ambas variables resultados <-by(bd, list(group1, group2), seqmeant)# Obtener todas las combinaciones posibles de los grupos combinaciones <-expand.grid(group1 =unique(group1), group2 =unique(group2), stringsAsFactors =FALSE)# Extraer los resultados y asociarlos con las combinaciones resultados_df <-do.call(rbind, lapply(seq_along(resultados), function(i) { group_name1 <-attr(resultados, "dimnames")[[1]][i] group_name2 <-attr(resultados, "dimnames")[[2]][i]data.frame(factor_inclusivo_1 = group_name1, factor_inclusivo_2 = group_name2, Mean = resultados[[i]]) }))# Unir los resultados con las combinaciones para rellenar los valores faltantes final_df <-merge(combinaciones, resultados_df, by.x =c("group1", "group2"), by.y =c("factor_inclusivo_1", "factor_inclusivo_2"), all.x =TRUE)return(final_df)}multinom_pivot_wider <-function(x) {# check inputs match expectatations# create tibble of results df <- tibble::tibble(outcome_level =unique(x$table_body$groupname_col)) df$tbl <- purrr::map( df$outcome_level,function(lvl) { gtsummary::modify_table_body( x, ~dplyr::filter(.x, .data$groupname_col %in% lvl) |> dplyr::ungroup() |> dplyr::select(-.data$groupname_col) ) } )tbl_merge(df$tbl, tab_spanner =paste0("**", df$outcome_level, "**"))}best_subset_multinom <-function(y, x.vars, data) {# y Nombre de la variable dependiente (cadena de texto)# x.vars Vector de nombres de predictores (caracter)# data Dataframe con los datos de entrenamiento# Cargar las librerías necesariasrequire(dplyr)require(purrr)require(tidyr)require(nnet)require(MASS)# Generar todas las combinaciones posibles de predictores predictors_list <-lapply(1:length(x.vars), function(i) {combn(x.vars, i, simplify =FALSE) }) |>unlist(recursive =FALSE)# Inicializar una lista para almacenar los resultados results <-list()# Iterar sobre cada combinación de predictoresfor (i inseq_along(predictors_list)) { predictors <- predictors_list[[i]] formula <-as.formula(paste(y, "~", paste(predictors, collapse ="+")))# Ajustar el modelo multinomial model <-tryCatch( nnet::multinom(formula, data = data, trace =FALSE),error =function(e) NULL )# Si el modelo se ajustó correctamente, almacenar los resultadosif (!is.null(model)) {# Extraer el AIC del modelo aic <-AIC(model)# Almacenar la información en una lista results[[length(results) +1]] <-list(predictors = predictors,model = model,AIC = aic ) } }# Convertir la lista de resultados en un dataframe results_df <- results |> purrr::map_df(function(res) {data.frame(predictors =paste(res$predictors, collapse ="+"),AIC = res$AIC,stringsAsFactors =FALSE ) })# Ordenar los modelos por AIC de menor a mayor results_df <- results_df |>arrange(AIC)return(results_df)}best_subset_multinom_interactions <-function(y, x.vars, data) {# y Nombre de la variable dependiente (cadena de texto)# x.vars Vector de nombres de predictores (caracter)# data Dataframe con los datos de entrenamiento# Cargar las librerías necesariasrequire(dplyr)require(purrr)require(tidyr)require(nnet)require(MASS)# Generar todas las combinaciones posibles de predictores (efectos principales) main_effects_list <-lapply(1:length(x.vars), function(i) {combn(x.vars, i, simplify =FALSE) }) |>unlist(recursive =FALSE)# Inicializar una lista para almacenar los resultados results <-list()# Iterar sobre cada combinación de efectos principalesfor (main_effects in main_effects_list) {# Generar términos de interacción de hasta 3 variables interaction_terms <-list()# Para interacciones de 2 variablesif (length(main_effects) >=2) { interaction_terms_2way <-combn(main_effects, 2, function(x) paste(x, collapse =":")) interaction_terms <-c(interaction_terms, interaction_terms_2way) }# Para interacciones de 3 variablesif (length(main_effects) >=3) { interaction_terms_3way <-combn(main_effects, 3, function(x) paste(x, collapse =":")) interaction_terms <-c(interaction_terms, interaction_terms_3way) }# Combinar efectos principales e interacciones all_terms <-c(main_effects, interaction_terms)# Generar todas las combinaciones posibles de términos (incluyendo interacciones)# Solo se incluyen interacciones si sus efectos principales están presentes term_combinations <-list()# Obtener todos los subconjuntos de efectos principales main_effects_subsets <-lapply(1:length(main_effects), function(i) {combn(main_effects, i, simplify =FALSE) }) |>unlist(recursive =FALSE)# Para cada subconjunto de efectos principalesfor (me in main_effects_subsets) {# Iniciar con los efectos principales terms <- me# Incluir interacciones solo si todos sus efectos principales están incluidos possible_interactions <- interaction_terms[sapply(interaction_terms, function(x) { vars_in_interaction <-unlist(strsplit(x, ":"))all(vars_in_interaction %in% me) }) ]# Generar todas las combinaciones de interacciones para incluir interaction_subsets <-list(NULL)if (length(possible_interactions) >0) { interaction_subsets <-lapply(1:length(possible_interactions), function(i) {combn(possible_interactions, i, simplify =FALSE) }) |>unlist(recursive =FALSE) }# Para cada combinación de interacciones, crear el conjunto completo de términosfor (ints in interaction_subsets) {if (is.null(ints)) { full_terms <- terms } else { full_terms <-c(terms, ints) }# Añadir a la lista de combinaciones de términos term_combinations <-append(term_combinations, list(full_terms)) } }# Ajustar modelos para cada combinación de términosfor (terms in term_combinations) { formula <-as.formula(paste(y, "~", paste(terms, collapse ="+")))# Ajustar el modelo multinomial model <-tryCatch( nnet::multinom(formula, data = data, trace =FALSE),error =function(e) NULL,warning =function(w) NULL )# Si el modelo se ajustó correctamente, almacenar los resultadosif (!is.null(model)) {# Extraer el BIC del modelo bic <-BIC(model)# Almacenar la información en la lista de resultados results[[length(results) +1]] <-list(predictors =paste(terms, collapse =" + "),model = model,BIC = bic ) } } }# Convertir la lista de resultados en un dataframe results_df <- results |> purrr::map_df(function(res) {data.frame(predictors = res$predictors,BIC = res$BIC,stringsAsFactors =FALSE ) })# Ordenar los modelos por BIC de menor a mayor results_df <- results_df |>arrange(BIC)return(results_df)}best_subset_multinom_interactions_parallel <-function(y, x.vars, data) {# y Nombre de la variable dependiente (cadena de texto)# x.vars Vector de nombres de predictores (caracter)# data Dataframe con los datos de entrenamiento# Cargar las librerías necesarias dentro de la funciónrequire(dplyr)require(purrr)require(tidyr)require(nnet)require(MASS)require(foreach)require(doParallel)require(progressr)# Iniciar los gestores de progresohandlers(global =TRUE)handlers("txt")# Generar todas las combinaciones posibles de predictores (efectos principales) main_effects_list <-lapply(1:length(x.vars), function(i) {combn(x.vars, i, simplify =FALSE) }) |>unlist(recursive =FALSE)# Inicializar una lista para almacenar las fórmulas de los modelos formulas_list <-list()# Generar todas las fórmulas posibles con interacciones hasta de 3 variablesfor (main_effects in main_effects_list) {# Generar términos de interacción de hasta 3 variables interaction_terms <-character(0) # Aseguramos que es un vector de caracteres# Para interacciones de 2 variablesif (length(main_effects) >=2) { interaction_terms_2way <-combn(main_effects, 2, function(x) paste(x, collapse =":"), simplify =TRUE) interaction_terms <-c(interaction_terms, interaction_terms_2way) }# Para interacciones de 3 variablesif (length(main_effects) >=3) { interaction_terms_3way <-combn(main_effects, 3, function(x) paste(x, collapse =":"), simplify =TRUE) interaction_terms <-c(interaction_terms, interaction_terms_3way) }# Generar todas las combinaciones posibles de efectos principales main_effects_subsets <-lapply(1:length(main_effects), function(i) {combn(main_effects, i, simplify =FALSE) }) |>unlist(recursive =FALSE)# Para cada subconjunto de efectos principalesfor (me in main_effects_subsets) {# Iniciar con los efectos principales terms <- me# Identificar interacciones cuyos efectos principales están en 'me'if (length(interaction_terms) >0) { possible_interactions <- interaction_terms[vapply(interaction_terms, function(x) { vars_in_interaction <-unlist(strsplit(x, ":"))all(vars_in_interaction %in% me) }, FUN.VALUE =logical(1)) ] } else { possible_interactions <-character(0) }# Generar todas las combinaciones posibles de estas interacciones interaction_subsets <-list(character(0)) # Incluir el caso sin interaccionesif (length(possible_interactions) >0) { interaction_combinations <-lapply(1:length(possible_interactions), function(i) {combn(possible_interactions, i, simplify =FALSE) }) |>unlist(recursive =FALSE) interaction_subsets <-c(interaction_subsets, interaction_combinations) }# Para cada combinación de interaccionesfor (ints in interaction_subsets) { full_terms <-c(terms, ints)# Crear la fórmula del modelo y almacenarla formula_str <-paste(y, "~", paste(full_terms, collapse ="+")) formulas_list <-append(formulas_list, list(formula_str)) } } }# Eliminar posibles duplicados de fórmulas formulas_list <-unique(formulas_list)# Total de modelos a ajustar total_models <-length(formulas_list)# Iniciar el progreso p <-progressor(steps = total_models)# Ajustar los modelos en paralelo usando foreach results_list <-foreach(i =1:total_models, .packages =c("nnet", "MASS"), .combine ='rbind') %dopar% { formula_str <- formulas_list[[i]] formula <-as.formula(formula_str)# Ajustar el modelo model <-tryCatch( nnet::multinom(formula, data = data, trace =FALSE),error =function(e) NULL,warning =function(w) NULL )# Actualizar el progresop(sprintf("Ajustando modelo %d de %d", i, total_models))# Si el modelo se ajustó correctamente, almacenar los resultadosif (!is.null(model)) { bic <-BIC(model)data.frame(predictors = formula_str,BIC = bic,stringsAsFactors =FALSE ) } else {NULL } }# Convertir los resultados a dataframe y ordenar por BIC results_df <-as.data.frame(results_list) results_df <- results_df |>arrange(BIC)return(results_df)}num_cores <- parallel::detectCores() -1cl <-makeCluster(num_cores)registerDoParallel(cl)#pacman job kableExtra tidyverse cluster WeightedCluster devtools TraMineR TraMineRextras NbClust haven ggseqplot gridExtra Tmisc factoextra reticulate withr rmarkdown quartooptions(knitr.kable.NA ='')#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:##:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#knitr::knit_hooks$set(time_it =local({ now <-NULLfunction(before, options) {if (before) {# record the current time before each chunk now <<-Sys.time() } else {# calculate the time difference after a chunk res <-ifelse(difftime(Sys.time(), now)>(60^2),difftime(Sys.time(), now)/(60^2),difftime(Sys.time(), now)/(60^1))# return a character string to show the time x<-ifelse(difftime(Sys.time(), now)>(60^2),paste("Tiempo que demora esta sección:", round(res,1), "horas"),paste("Tiempo que demora esta sección:", round(res,1), "minutos"))paste('<div class="message">', gsub('##', '\n', x),'</div>', sep ='\n') } }}))knitr::opts_chunk$set(time_it =TRUE)#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:format_cells <-function(df, rows ,cols, value =c("italics", "bold", "strikethrough")){# select the correct markup# one * for italics, two ** for bold map <-setNames(c("*", "**", "~~"), c("italics", "bold", "strikethrough")) markup <- map[value] for (r in rows){for(c in cols){# Make sure values are not factors df[[c]] <-as.character( df[[c]])# Update formatting df[r, c] <-ifelse(nchar(df[r, c])==0,"",paste0(markup, gsub(" ", "", df[r, c]), markup)) } }return(df)}#To produce line breaks in messages and warningsknitr::knit_hooks$set(error =function(x, options) {paste('\n\n<div class="alert alert-danger" style="font-size: small !important;">',gsub('##', '\n', gsub('^##\ Error', '**Error**', x)),'</div>', sep ='\n') },warning =function(x, options) {paste('\n\n<div class="alert alert-warning" style="font-size: small !important;">',gsub('##', '\n', gsub('^##\ Warning:', '**Warning**', x)),'</div>', sep ='\n') },message =function(x, options) {paste('<div class="message" style="font-size: small !important;">',gsub('##', '\n', x),'</div>', sep ='\n') })#_#_#_#_#_#_#_#_#_#_#_#_#_invisible("Function to format CreateTableOne into a database")as.data.frame.TableOne <-function(x, ...) {capture.output(print(x,showAllLevels =TRUE, varLabels = T,...) -> x) y <-as.data.frame(x) y$characteristic <- dplyr::na_if(rownames(x), "") y <- y |>fill(characteristic, .direction ="down") |> dplyr::select(characteristic, everything())rownames(y) <-NULL y}#_#_#_#_#_#_#_#_#_#_#_#_#_# Austin, P. C. (2009). The Relative Ability of Different Propensity # Score Methods to Balance Measured Covariates Between # Treated and Untreated Subjects in Observational Studies. Medical # Decision Making. https://doi.org/10.1177/0272989X09341755smd_bin <-function(x,y){ z <- x*(1-x) t <- y*(1-y) k <-sum(z,t) l <- k/2return((x-y)/sqrt(l))}#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:format_table_vec <-function(tbl, digits =1) { counts <-as.numeric(tbl) percentages <-prop.table(tbl) *100 formatted <-sapply(seq_along(counts), function(i) { p_val <- percentages[i]# Si el porcentaje es prácticamente entero, formatea sin decimalesif (abs(p_val -round(p_val)) < .Machine$double.eps^0.5) { p_str <-sprintf("%.0f", p_val) } else { p_str <-sprintf(paste0("%.", digits, "f"), p_val) }paste0(counts[i], " (", p_str, ")") }) formatted}#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:if(.Platform$OS.type =="windows") withAutoprint({memory.size()memory.size(TRUE)memory.limit()})
> memory.size()
[1] Inf
> memory.size(TRUE)
[1] Inf
> memory.limit()
[1] Inf
Código
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:func_tab_range_clus<-function(range_clus){rbind.data.frame(lapply(list(as.vector(rev(sort(table(range_clus$clustering$cluster2)))),as.vector(rev(sort(table(range_clus$clustering$cluster3)))),as.vector(rev(sort(table(range_clus$clustering$cluster4)))),as.vector(rev(sort(table(range_clus$clustering$cluster5)))),as.vector(rev(sort(table(range_clus$clustering$cluster6)))),as.vector(rev(sort(table(range_clus$clustering$cluster7)))),as.vector(rev(sort(table(range_clus$clustering$cluster8)))),as.vector(rev(sort(table(range_clus$clustering$cluster9)))),as.vector(rev(sort(table(range_clus$clustering$cluster10)))),as.vector(rev(sort(table(range_clus$clustering$cluster11)))),as.vector(rev(sort(table(range_clus$clustering$cluster12)))),as.vector(rev(sort(table(range_clus$clustering$cluster13)))),as.vector(rev(sort(table(range_clus$clustering$cluster14)))),as.vector(rev(sort(table(range_clus$clustering$cluster15)))) ),function(x) { length_out <-max(sapply(list(as.vector(rev(sort(table(range_clus$clustering$cluster2)))),as.vector(rev(sort(table(range_clus$clustering$cluster3)))),as.vector(rev(sort(table(range_clus$clustering$cluster4)))),as.vector(rev(sort(table(range_clus$clustering$cluster5)))),as.vector(rev(sort(table(range_clus$clustering$cluster6)))),as.vector(rev(sort(table(range_clus$clustering$cluster7)))),as.vector(rev(sort(table(range_clus$clustering$cluster8)))),as.vector(rev(sort(table(range_clus$clustering$cluster9)))),as.vector(rev(sort(table(range_clus$clustering$cluster10)))),as.vector(rev(sort(table(range_clus$clustering$cluster11)))),as.vector(rev(sort(table(range_clus$clustering$cluster12)))),as.vector(rev(sort(table(range_clus$clustering$cluster13)))),as.vector(rev(sort(table(range_clus$clustering$cluster14)))),as.vector(rev(sort(table(range_clus$clustering$cluster15)))) ), length))c(x, rep(NA, length_out -length(x))) } ))|>t() |>data.frame()|>`rownames<-`(NULL)}frobenius_norm <-function(matrix1, matrix2) {if (!all(dim(matrix1) ==dim(matrix2))) {stop("Matrices must have the same dimensions") }# Replace NA values with 0 (or any other desired default) matrix1[is.na(matrix1)] <-0 matrix2[is.na(matrix2)] <-0# Calculate the residuals residuals <- matrix1 - matrix2# Frobenius norm frobenius <-sqrt(sum(residuals^2))return(frobenius)}#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:confcqi2 <-function(nullstat, quant, n){ alpha <- (1-quant)/2#calpha <- alpha+(alpha-1)/n#print(c(calpha, alpha))#minmax <- quantile(nullstat, c(calpha, 1-calpha)) minmax <-quantile(nullstat, c(alpha, 1-alpha))return(minmax)}normstatcqi2 <-function(bcq, stat, norm=TRUE){ origstat <- bcq$clustrange$stats[, stat] nullstat <- bcq$stats[[stat]]#normstat <- rbind(nullstat, origstat)if(norm){for(i inseq_along(origstat)){ mx <-mean(nullstat[, i]) sdx <-sd(nullstat[, i]) nullstat[ , i] <- (nullstat[, i]-mx)/sdx origstat[i] <- (origstat[i]-mx)/sdx } } alldatamax <-apply(nullstat, 1, max)#as.vector(xx) sumcqi <-list(origstat=origstat, nullstat=nullstat, alldatamax=alldatamax)return(sumcqi)}print.seqnullcqi.powder <-function(x, norm =FALSE, quant =0.95, digits =2, append =FALSE, ...) {cat("Parametric bootstrap cluster analysis validation\n")cat("Sequence analysis null model:", deparse(x$nullmodel), "\n")cat("Number of bootstraps:", x$R, "\n")cat("Clustering method:", ifelse(x$kmedoid, "PAM/K-Medoid", paste0("hclust with ", x$hclust.method)), "\n")cat("Seqdist arguments:", deparse(x$seqdist.args), "\n\n\n") alls <-as.data.frame(x$clustrange$stats) quants <-rep("", ncol(alls))names(quants) <-colnames(alls)for (ss incolnames(alls)) { sumcqi <-normstatcqi2(x, stat = ss, norm = norm) alls[, ss] <-as.character(round(sumcqi$origstat, digits = digits)) borne <-as.character(round(confcqi2(sumcqi$alldatamax, quant, x$R), digits = digits)) quants[ss] <-paste0("[", borne[1], "; ", borne[2], "]") } results_tibble <- tibble::as_tibble(rbind(alls, rep("", length(quants)), quants))# Print a summary to the console for immediate feedbackrownames(results_tibble) <-c(rownames(x$clustrange$stats), "", paste("Null Max-T", quant, "interval")) results_df <-as.data.frame(results_tibble)print(results_tibble, ...)return(list(results_tibble= results_tibble, results_df= results_df ))}#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:##:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:# Función para aplicar la prueba de Fisher a todas las combinaciones de filas usando todas las columnasfisher_posthoc_all_cols <-function(contingency_table) {# Obtener combinaciones de filas (pares) row_pairs <-combn(rownames(contingency_table), 2, simplify =FALSE)# Aplicar la prueba de Fisher a cada par de filas usando todas las columnas al mismo tiempo results <-map_dfr(row_pairs, function(pair) {# Crear tabla de 2xN para el par de filas en todas las columnas sub_table <- contingency_table[pair, , drop =FALSE]# Aplicar el test de Fisher test_result <-fisher.test(sub_table, simulate.p.value=T,B=1e4)# Devolver los resultados en un data frametibble(Row1 = pair[1],Row2 = pair[2],p.value = test_result$p.value ) })# Ajustar p-valores usando el método de Holm results <- results |>mutate(p.adjusted =p.adjust(p.value, method ="holm"))return(results)}#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:save_base_plot_as_grob <-function(plot_expr, res=300, width =1600, height=1200) {# Crea un archivo temporal con extensión .png filename <-tempfile(fileext =".png")# Guarda el gráfico en alta resolución en el archivo temporalpng(filename, width = width, height = height, res = res)replayPlot(plot_expr) # Reproduce el gráfico grabadodev.off() # Cierra el dispositivo gráfico# Convierte el archivo PNG en un objeto gráfico (grob) grob <- grid::rasterGrob(png::readPNG(filename), interpolate =TRUE)return(grob) # Devuelve el grob}#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:chisq_cramerv<-function(contingency_table){ chisq_test <-chisq.test(contingency_table) cramers_v <-sqrt(chisq_test$statistic / (sum(contingency_table) * (min(dim(contingency_table)) -1)))list(chisq_statistic=sprintf("%1.2f", chisq_test$statistic), chisq_df= chisq_test$parameter, chisq_p_value =ifelse(chisq_test$p.value<.001, "<0.001", sprintf("%1.4f", chisq_test$p.value)), cramers_v =sprintf("%1.2f", cramers_v))}#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#oneway_anova_effect_size <-function(values, group) {# Perform one-way ANOVA anova_result <-aov(values ~ group)# Summarize ANOVA results anova_summary <-summary(anova_result)# Extract sums of squares ss_between <- anova_summary[[1]]$"Sum Sq"[1] ss_total <-sum(anova_summary[[1]]$"Sum Sq")# Calculate eta-squared eta_squared <- ss_between / ss_total# Return ANOVA summary and effect sizelist(anova_summary = anova_summary,eta_squared = eta_squared )}#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#duración en cada estadoseq_mean_t<-function(bd=NULL, group= group){ resultados<-by(bd, group, seqmeant)do.call(rbind, lapply(names(resultados), function(name) {data.frame(factor_inclusivo = name, resultados[[name]]) }))}seqtrate_t<-function(bd=NULL, group= group){# Utilizar la función 'by' para calcular las tasas agrupadas por 'glosa_sexo' resultados <-by(bd, group, seqtrate)# Convertir los resultados en un data frame en formato largo resultados_long <-do.call(rbind, lapply(names(resultados), function(sexo) { df <-as.data.frame(resultados[[sexo]]) df$from <-rownames(df) df$glosa_sexo <- sexo df }))# Usar tidyr para convertir a formato largolibrary(tidyr) resultados_long <-pivot_longer(resultados_long, cols =-c(from, glosa_sexo), names_to ="to", values_to ="rate")# Mostrar el data frame finalprint(resultados_long)}seqcount_t<-function(bd=NULL, group= group){# Utilizar la función 'by' para calcular las tasas agrupadas por 'glosa_sexo' resultados <-by(bd, group, function(x) seqtrate(x, count =TRUE))# Convertir los resultados en un data frame en formato largo resultados_long <-do.call(rbind, lapply(names(resultados), function(sexo) { df <-as.data.frame(resultados[[sexo]]) df$from <-rownames(df) df$glosa_sexo <- sexo df }))# Usar tidyr para convertir a formato largolibrary(tidyr) resultados_long <-pivot_longer(resultados_long, cols =-c(from, glosa_sexo), names_to ="to", values_to ="count")# Mostrar el data frame finalprint(resultados_long)}#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
Descartamos objetos que no serán de interés para el análisis
Código
# Define el vector con los nombres de los objetos que quieres conservarmis_objetos <-c("as.data.frame.TableOne", "avs_por_cluster_month", "avs_por_cluster_quarter","best_subset_multinom", "best_subset_multinom_interactions", "best_subset_multinom_interactions_parallel", "categories_hac_om2_m", "categories_pam_om4_q", "chisq_cramerv", "cl", "costmatrix_month", "costmatrix_quarter", "data_long_establecimiento_2024_std", "df_filled2", "dist_month_lcs", "dist_month_om", "dist_quarter_lcs", "dist_quarter_om", "dt_ing_calendar_month_t_desde_primera_adm_dedup", "dt_ing_calendar_quarter_t_desde_primera_adm_dedup", "fisher_posthoc_all_cols", "folder_path", "format_cells", "format_table_vec", "ing_dt_ing_calendar_month_t_desde_primera_adm_dedup_wide2_cens", "ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens", "lcs_dist_month", "lcs_dist_month_c", "lcs_dist_quarter", "lcs_dist_quarter_c", "max_clusters", "multinom_pivot_wider", "new_labels", "new_labels2", "num_cores", "om_dist_month", "om_dist_month_c", "om_dist_quarter", "om_dist_quarter_c", "oneway_anova_effect_size", "pamRange_month_lcs", "pamRange_month_lcs2", "pamRange_month_om", "pamRange_month_om2", "pamRange_quarter_lcs", "pamRange_quarter_lcs2", "pamRange_quarter_om", "pamRange_quarter_om2", "ratio_plot", "resultados_list", "seq_plot_hc_om2_m", "seq_plot_hc_om3_m", "seq_plot_pam_om2_m", "seq_plot_pam_om2_q", "seq_plot_pam_om3_q", "seq_plot_pam_om4_q", "seq_plot_pam_om6_q", "sil_hc_om_clus2_m", "sil_hc_om_clus2_q", "sil_hc_om_clus3_m", "sil_pam_om_clus2_m", "sil_pam_om_clus2_q", "sil_pam_om_clus3_q", "sil_pam_om_clus4_q", "sil_pam_om_clus6_q", "smd_bin", "States_Wide.seq_month_t_prim_adm", "States_Wide.seq_month_t_prim_adm_cens", "States_Wide.seq_month_t_prim_adm_RM_cens", "States_Wide.seq_quarter_t_prim_adm", "States_Wide.seq_quarter_t_prim_adm_cens", "States_Wide.seq_quarter_t_prim_adm_RM_cens", "seqcount_t", "seq_mean_t", "seqtrate_t")# Filtrar el entorno global y eliminar todo lo que NO esté en ese vectorrm(list =setdiff(ls(), mis_objetos)); gc()
used (Mb) gc trigger (Mb) max used (Mb)
Ncells 4482668 239.5 16486446 880.5 20608057 1100.6
Vcells 358249052 2733.3 2136999950 16304.1 2569572024 19604.3
Tiempo que demora esta sección: 0 minutos
Resultados
1. Trimestral
1.1. PAM (OM), sol 4 cluster- diagnósticos
Código
invisible("Me da buena: 0,61 en promedio. El problema está con 5710 y 6036 que son pésimos")sil_pam_om_clus4_q_nostd<-silhouette(as.integer(pamRange_quarter_om$clustering$cluster4), as.dist(dist_quarter_om))# Crear etiquetas personalizadascluster_labels <-paste0("Cluster ", seq_along(attr(summary(sil_pam_om_clus4_q_nostd)$clus.avg.widths, "dimnames")[[1]]), ":\nAWS ", sprintf("%1.2f",summary(sil_pam_om_clus4_q_nostd)$clus.avg.widths))# Graficar con etiquetas personalizadasfviz_silhouette( sil_pam_om_clus4_q_nostd, lab.clusters = cluster_labels, # Etiquetas personalizadas para los clústeresprint.summary=F) +scale_fill_grey(start =0.2, end =0.8, labels = cluster_labels) +# Escala de grisesscale_color_grey(start =0.2, end =0.8, labels = cluster_labels)+# Escala de grises para los bordeaggtitle(NULL)+labs(y="Ancho medio de la silueta", x="Conglomerados")# Elimina el título
cat("de quienes pertenecen al grupo de un semestre")
de quienes pertenecen al grupo de un semestre
Código
length(subset(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens, clus_pam_om4=="6522, Un semestre TSM(1)")$run)
[1] 399
Código
diag_pam_om4_6522<-df_filled2 |> dplyr::filter(run %in%subset(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens, clus_pam_om4=="6522, Un semestre TSM(1)")$run) |> dplyr::select(run, diag1, diag2, diag3, diag4, diag5, diag6, diag7, diag8, diag9, diag10, diag11, fecha_egreso_rec_fmt, estab_homo) |> dplyr::group_by(run) |> dplyr::filter(row_number() !=1) |># Elimina la primera observación de cada run dplyr::mutate(all_diags =paste(na.omit(c(diag1, diag2, diag3, diag4, diag5, diag6, diag7, diag8, diag9, diag10, diag11)), collapse =", ") ) |> dplyr::summarise(all_diags =first(all_diags),fecha_egreso_rec_fmt =first(fecha_egreso_rec_fmt),estab_homo =first(estab_homo) ) |> dplyr::ungroup() |> dplyr::pull(all_diags) |># Extraer la columna all_diags como vectorstrsplit(split =", ") |># Separar cada diagnóstico por comasunlist() # Convertir la lista en un vectorinvisible("head(arrange(data.frame(table(diag_pam_om4_6522)),-Freq),20)")invisible("Para chatgpt= estos son códigos de CIE-10, descríbeme brevemente cada uno en markdown en formato 'Cód. CIE-10 (n=Freq) - [descripción] '")cat("número de diagnósticos distintos")
## F329 (n=134) - Episodio depresivo no especificado # F322 (n=103) - Episodio depresivo grave sin síntomas psicóticos # F609 (n=96) - Otros trastornos específicos de la personalidad # F603 (n=95) - Trastorno de la personalidad emocionalmente inestable, tipo límite # F209 (n=66) - Esquizofrenia no especificada # F192 (n=53) - Trastornos mentales y del comportamiento debidos al uso de múltiples drogas y al uso de otras sustancias psicoactivas # F319 (n=49) - Trastorno depresivo recurrente, episodio actual no especificado # F200 (n=45) - Esquizofrenia paranoide # F432 (n=38) - Trastornos de adaptación # F323 (n=35) - Episodio depresivo grave con síntomas psicóticos # Z915 (n=32) - Antecedentes personales de traumatismo no clasificado en otra parte # C490 (n=29) - Neoplasia maligna del tejido conectivo y de los tejidos blandos, de localización no especificada # G909 (n=29) - Trastorno del sistema nervioso no especificado # F431 (n=28) - Reacción al estrés agudo y trastornos de adaptación # E669 (n=26) - Obesidad, no especificada# Z511 (n=24) - Atención sanitaria para radioterapia # F29X (n=21) - Psicosis no orgánica no especificada # F419 (n=20) - Trastorno de ansiedad no especificado # C901 (n=16) - Leucemia de células plasmáticas# F199 (n=16) - Trastornos mentales y de compurtamiento por el uso de múltiples drogas y uso de otras sustancias piscoactivas
Tiempo que demora esta sección: 0 minutos
Entre quienes se encuentren en un semestre en el sistema por TSM y presentan un segundo episodio, de las 20 condiciones más frecuentes, al menos el 50% se caracteriza por episodios depresivos no especificado (F329), trastornos de la personalidad tipo límite (F603), otros trastornos específicos de la personalidad (F609) y episodios depresivos graves sin síntomas psicóticos (F322) y esquizofrenia no especificada (F209).
Código
# Trastornos mentales orgánicos: F00.0-F09.9 (F000 a F099)in_organic <-Vectorize(function(codigo) {if (is.na(codigo)) return(FALSE)grepl("^F0[0-9]", codigo)})# Trastornos por uso de sustancias: F10.0-F19.9 (F10 a F19)in_substance <-Vectorize(function(codigo) {if (is.na(codigo)) return(FALSE)grepl("^F1[0-9]", codigo)})# Esquizofrenia: F20.0-F20.9 (F20)in_esquizofrenia <-Vectorize(function(codigo) {if (is.na(codigo)) return(FALSE)grepl("^F20", codigo)})# Otros trastornos psicóticos no afectivos: F21.0-F29.9 (F21 a F29)in_otros_psicoticos <-Vectorize(function(codigo) {if (is.na(codigo)) return(FALSE)grepl("^F2[1-9]", codigo)})# Trastornos bipolares: F30.0-F31.9 (F30 o F31)in_bipolares <-Vectorize(function(codigo) {if (is.na(codigo)) return(FALSE)grepl("^(F30|F31)", codigo)})# Trastornos depresivos y otros del estado de ánimo: F32.0-F39.9 (F32 a F39)in_depresivos <-Vectorize(function(codigo) {if (is.na(codigo)) return(FALSE)grepl("^F3[2-9]", codigo)})# Trastornos de ansiedad: F40.0-F49.9 (F40 a F49)in_ansiedad <-Vectorize(function(codigo) {if (is.na(codigo)) return(FALSE)grepl("^F4[0-9]", codigo)})# Trastornos de la personalidad: F60.0-F69.9 (F60 a F69)in_personalidad <-Vectorize(function(codigo) {if (is.na(codigo)) return(FALSE)grepl("^F6[0-9]", codigo)})# LESIONES AUTOINFLIGIDAS INTENCIONALMENTE (X60-X84)in_lesiones <-Vectorize(function(codigo) {if (is.na(codigo)) return(FALSE)grepl("^X6[0-9]|^X7[0-9]|^X8[0-4]", codigo)})#Mortier, P., Conde, S., Alayo, I., Amigo, F., Ballester, L., Amell, R. C., Guinart, D., Contaldo, S. F., Ferrer, M., Leis, A., Mayer, M. A., Portillo-Van Diest, A., Puértolas-Gracia, B., Ramírez-Anguita, J. M., Peña-Salazar, C., Sanz, F., Kessler, R. C., Palao, D., Sola, V. P., . . . Alonso, J. (2024). Premature Death, Suicide, and Nonlethal Intentional Self-Harm After Psychiatric Discharge. JAMA Network Open, 7(6), e2417131. https://doi.org/10.1001/jamanetworkopen.2024.17131in_lesiones_amplio <-Vectorize(function(codigo) {if (is.na(codigo)) return(FALSE)grepl("^(T3[6-9]|T4[0-9]|T5[0-9]|T6[0-5])", codigo)})agregar_columnas_icd <-function(df) {# Especificamos las columnas de diagnóstico diag_cols <-paste0("diag", 1:11) df |>mutate(organic =rowSums(across(all_of(diag_cols), ~in_organic(.))) >0,substance =rowSums(across(all_of(diag_cols), ~in_substance(.))) >0,esquizofrenia =rowSums(across(all_of(diag_cols), ~in_esquizofrenia(.))) >0,otros_psicoticos =rowSums(across(all_of(diag_cols), ~in_otros_psicoticos(.))) >0,bipolares =rowSums(across(all_of(diag_cols), ~in_bipolares(.))) >0,depresivos =rowSums(across(all_of(diag_cols), ~in_depresivos(.))) >0,ansiedad =rowSums(across(all_of(diag_cols), ~in_ansiedad(.))) >0,personalidad =rowSums(across(all_of(diag_cols), ~in_personalidad(.))) >0 )}#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_cat("2025-03-01: general")
#[1] 1687190invisible("Quienes tienen un registro RSH o en el histórico una autoidentificación por MINSAL")ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$ppoo_minsal_y_rsh_2010 <-ifelse(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$run %in% filtered_df_csv$run,1,0)rm(filtered_df_csv)ci_cya_conadi<-readr::read_csv("cya_conadi_o_ci_conadi.csv")nrow(ci_cya_conadi)
invisible("SI tiene un registro CONADI CyA o CI, entonces es 1")ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$ci_cya_conadi <-ifelse(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$run %in% ci_cya_conadi$RUN,1,0)cat("Tabla de interesección entre pueblos originarios y CONADI, histórico (2005-)\n")
Tabla de interesección entre pueblos originarios y CONADI, histórico (2005-)
print("Estaríamos elimnando un 18% de casos de registros en otro momento")
[1] "Estaríamos elimnando un 18% de casos de registros en otro momento"
Código
glosa_pueblo_originario_rec_num<-as.numeric(factor(glosa_pueblo_originario_rec))-1seqcompare_ppoo_quarter_om_ppoo_rec<-TraMineRextras::seqCompare(States_Wide.seq_quarter_t_prim_adm, group=glosa_pueblo_originario_rec_num, method="OM", sm=costmatrix_quarter, seed=2125, s=1e1, stat="all")# s= Integer. Default 100. The size of random samples of sequences. When 0, no sampling is done.
Tiempo que demora esta sección: 0.2 minutos
Generamos un gráfico de PPOO por cada conglomerado.
Código
ppoo_clus_pre_pam_om4_q<- df_filled2[,c("run","glosa_pueblo_originario")] |> dplyr::left_join(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens[,c("run", "clus_pam_om4","ci_cya_conadi", "ppoo_minsal_y_rsh_2010")], by="run", multiple="first") |> dplyr::mutate(glosa_pueblo_originario_rec= dplyr::case_when(glosa_pueblo_originario =="NINGUNO"& (ci_cya_conadi ==1| ppoo_minsal_y_rsh_2010==1)~"DESCONOCIDO", T~glosa_pueblo_originario)) |> janitor::tabyl(glosa_pueblo_originario_rec, clus_pam_om4) |> janitor::adorn_percentages("row")invisible("2025-04-13")cat("Sólo RAPA NUI y MAPUCHE (ver por RUN, al menos una observación que se identifique RAPA NUI o MAPUCHE)\n")
Sólo RAPA NUI y MAPUCHE (ver por RUN, al menos una observación que se identifique RAPA NUI o MAPUCHE)
Fisher's Exact Test for Count Data with simulated p-value (based on
1e+05 replicates)
data: ppoo_clus_pre_pam_om4_q_clasificaciones
p-value = 0.0102
alternative hypothesis: two.sided
Fisher's Exact Test for Count Data with simulated p-value (based on
1e+05 replicates)
data: ppoo_clus_pre_pam_om4_q_clasificaciones_sin_asw_neg
p-value = 0.00237
alternative hypothesis: two.sided
Código
#p-value = 0.00258cat("No vale la pena repetir el post-hoc, es muy similar todo (lo hice)")
No vale la pena repetir el post-hoc, es muy similar todo (lo hice)
Código
#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#Lcat("el primer evento ver qué reportaron a MINSAL en PPOO")
el primer evento ver qué reportaron a MINSAL en PPOO
Fisher's Exact Test for Count Data with simulated p-value (based on
1e+05 replicates)
data: ppoo_clus_pre_pam_om4_q_clasificaciones2
p-value = 0.01382
alternative hypothesis: two.sided
cat("Tampoco resulta significativo en https://adaptivedesignstrial.shinyapps.io/posthoc/\n")
Tampoco resulta significativo en https://adaptivedesignstrial.shinyapps.io/posthoc/
Código
print("Eso si resulta interesante que los 25 usuarios desconcidos puro son 0.02 aunque en conjunto no hay bonferroni ni nada, pero son mas grandes; también es más grande Mapuche")
[1] "Eso si resulta interesante que los 25 usuarios desconcidos puro son 0.02 aunque en conjunto no hay bonferroni ni nada, pero son mas grandes; también es más grande Mapuche"
Código
#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#L#Lcat("Ahora descartando siluetas negativas, para el primer evento ver qué reportaron a MINSAL en PPOO")
Ahora descartando siluetas negativas, para el primer evento ver qué reportaron a MINSAL en PPOO
Fisher's Exact Test for Count Data with simulated p-value (based on
1e+05 replicates)
data: ppoo_clus_pre_pam_om4_q_clasificaciones_sin_asw_neg2
p-value = 0.00249
alternative hypothesis: two.sided
Código
#p-value = 0.00266cat("No vale la pena repetir el post-hoc, es muy similar todo (lo hice)")
No vale la pena repetir el post-hoc, es muy similar todo (lo hice)
Tiempo que demora esta sección: 0 minutos
Código
cat("(__________________________________________________________________)\n")reshape2::melt(ppoo_clus_pre_pam_om4_q, id.vars ="glosa_pueblo_originario_rec") |> dplyr::mutate(glosa_pueblo_originario_rec= dplyr::recode(glosa_pueblo_originario_rec, "OTRO (ESPECIFICAR)"="OTRO(n=86)", #"OTRO(n=77)", "RAPA NUI (PASCUENSE)"="RAPA NUI(n=37)", #"RAPA NUI(n=34)", "YAGÁN (YÁMANA)"="YAGÁN(n=2)",#"YAGÁN(n=2)","AYMARA"="AYMARA(n=15)",#"AYMARA(n=13)","COLLA"="COLLA(n=6)",#"COLLA(n=6)","DIAGUITA"="DIAGUITA(n=3)",#"DIAGUITA(n=3)","KAWÉSQAR"="KAWÉSQAR(n=4)",#"KAWÉSQAR(n=4)","MAPUCHE"="MAPUCHE(n=299)",#"MAPUCHE(n=255)","DESCONOCIDO"=".DESCONOCIDO(n=2.353)",#".DESCONOCIDO(n=1.985)","NINGUNO"=".NINGUNO(n=10.425)"#".NINGUNO(n=9.156)" )) |>ggplot(aes(x = glosa_pueblo_originario_rec, y = value, fill = variable)) +geom_bar(stat ="identity", position ="fill") +scale_fill_manual(values =c("#D2B48C", "#E27A5B", "#708090", "#6B8E23")) +labs(title =NULL,x ="Grupo Étnico",y ="Proporción de reportes",fill ="Grupos") +# Cambia el título de la leyenda a "Grupos"theme_minimal() +theme(axis.text.y =element_text(size =12), # Tamaño de las etiquetas de los grupos étnicosaxis.text.x =element_text(size =12), # Tamaño de las etiquetas del eje Xaxis.title.x =element_text(size =14), # Tamaño del título del eje Xaxis.title.y =element_text(size =14), # Tamaño del título del eje Yplot.title =NULL, # Tamaño y estilo del título del gráficolegend.title =element_text(size =14, margin =margin(b =-.1)), # Tamaño del título de la leyendalegend.spacing.y =unit(1.5, "lines"),legend.box.spacing =unit(0.5, "lines"), # Controla el espacio entre la leyenda y el gráficolegend.margin =margin(5, 5, 5, 5), legend.key.height =unit(1, "cm"), legend.text =element_text(size =12) # Tamaño del texto de la leyenda ) +coord_flip() # Hacer el gráfico horizontalggsave("_figs/grafico_ancho_achatado_pam_om4_q_25.png", width =10, height =5, dpi=1000)ppoo_clus_pre_pam_om4_q_rapanui<- df_filled2[,c("run","glosa_pueblo_originario")] |> dplyr::left_join(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens[,c("run", "clus_pam_om4","factor_inclusivo_real_hist_mas_autperc")], by="run", multiple="first") |> dplyr::mutate(glosa_pueblo_originario_rec= dplyr::case_when(glosa_pueblo_originario=="NINGUNO"& factor_inclusivo_real_hist_mas_autperc!="00"~"DESCONOCIDO", glosa_pueblo_originario=="NINGUNO"~"NINGUNO", glosa_pueblo_originario=="RAPA NUI (PASCUENSE)"~"RAPA NUI", T~"RESTO")) |> janitor::tabyl(glosa_pueblo_originario_rec, clus_pam_om4) |> janitor::adorn_percentages("row")cat("\norigen, en caso que sospechemos de los patrones de hospitalización observados para RAPA NUI \n")df_filled2[,c("run","glosa_pueblo_originario")] |> dplyr::left_join(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens[,c("run", "clus_pam_om4","factor_inclusivo_real_hist_mas_autperc", "codigo_region_rec_base")], by="run", multiple="first")|> dplyr::mutate(glosa_pueblo_originario_rec= dplyr::case_when(glosa_pueblo_originario=="NINGUNO"& factor_inclusivo_real_hist_mas_autperc!="00"~"DESCONOCIDO", glosa_pueblo_originario=="NINGUNO"~"NINGUNO", glosa_pueblo_originario=="RAPA NUI (PASCUENSE)"~"RN", T~"RESTO")) |> dplyr::filter(glosa_pueblo_originario_rec=="RN") |> janitor::tabyl(codigo_region_rec_base)
(__________________________________________________________________)
origen, en caso que sospechemos de los patrones de hospitalización observados para RAPA NUI
codigo_region_rec_base n percent
RM 22 0.5945946
noRM 15 0.4054054
Eventos hospitalarios y reportes de pertenencia a PPOO MINSAL, por conglomerados
Tiempo que demora esta sección: 0 minutos
1.1.1. Trayectorias
Vemos los gráficos de las trayectorias
Código
categories_pam_om4_q<-attr(States_Wide.seq_quarter_t_prim_adm_cens, "labels")new_labels <- categories_pam_om4_qnew_labels[which(categories_pam_om4_q =="Otras causas")] <-"Otras\ncausas"#new_labels[which(categories == "Consumo\nde sustancias")] <- "Consumo de\nsustancias"# Creamos un vector con las columnas llenando con NA si faltan valoressil_pam_om_clus4_q <-wcSilhouetteObs(as.dist(dist_quarter_om), pamRange_quarter_om$clustering$cluster4, measure="ASW")seq_plot_pam_om4_q <-ggseqiplot(States_Wide.seq_quarter_t_prim_adm_cens, group= ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$clus_pam_om4,facet_ncol=2, facet_nrow=2, sortv=sil_pam_om_clus4_q) +theme(legend.position ="none")+labs(x="Trimestres", y="# IDs de usuarios")+#guides(fill = guide_legend(nrow = 1))+theme(panel.spacing =unit(0.1, "lines"), # Reduce el espaciado entre los panelesaxis.text.y =element_text(size =15), # Tamaño de las etiquetas de los grupos étnicosaxis.text.x =element_text(size =15), # Tamaño de las etiquetas del eje Xaxis.title.x =element_text(size =15), # Tamaño del título del eje Xaxis.title.y =element_text(size =15, margin =margin(r =-10)),#,margin = margin(l = -10)),strip.text =element_text(size =11, margin =margin(b =-15)),legend.text =element_text(size =15),legend.spacing.x =unit(0.1, 'cm'), # Alinea el título de la leyenda hacia la izquierdalegend.box.margin =margin(t =0, r =0, b =0, l =-50),legend.position ="bottom", legend.justification ="left",panel.spacing.y =unit(0.5, "lines"),strip.placement ="outside", # Para colocar las tiras fuera de los ejesstrip.background =element_blank() # Elimina el fondo para que parezca más espacioso#legend.key.size = unit(1.5, "lines"), # Aumenta el tamaño de los símbolos en la leyenda )+guides(fill =guide_legend(nrow =1)) +scale_fill_manual(labels = new_labels, values=c("#E2725B", "#556B2F", "#D2B48C",#"#8B4513","#FFFFFF","#808080","#000000"))+scale_color_manual(labels = new_labels, values=c("#E2725B", "#556B2F", "#D2B48C",#"#8B4513","#FFFFFF","#808080","#000000"))seq_plot_pam_om4_q ggsave(filename="_figs/clusters_pam_om4_q_mod_25.png", seq_plot_pam_om4_q, width =9.5, height =5.5, dpi=1000)#:#:#:#:#:#:#:#:#:#:#:ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$clus_pam_om4_eng <- ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$clus_pam_om4levels(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$clus_pam_om4_eng) <-c("1.First quarter, MH(6623)", "2.First quarter, SUD(6612)", "3.First semester, MH(6522)", "4.First quarter, comorbidity(6574)")new_labels_eng <-c("Substance\nuse disorder","Comorbidity","Mental health\ndisorder","Not in\n hospital","Other\ncauses","Censored")# Plot with translated labels and Times New Roman fontseq_plot_pam_om4_q_eng <-ggseqiplot(States_Wide.seq_quarter_t_prim_adm_cens,group = ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$clus_pam_om4_eng,facet_ncol =2, facet_nrow =2, sortv = sil_pam_om_clus4_q) +labs(x ="Quarters Since First MH/SUD Hospitalization", y ="Patient index number") +theme(text =element_text(family ="Times"),panel.spacing =unit(0.1, "lines"),axis.text.y =element_text(size =15),axis.text.x =element_text(size =15),axis.title.x =element_text(size =15),axis.title.y =element_text(size =15, margin =margin(r =-8)),strip.text =element_text(size =15, margin =margin(b =-15)),legend.text =element_text(size =15),legend.spacing.x =unit(0.1, 'cm'),legend.box.margin =margin(t =0, r =0, b =0, l =-50),legend.position ="bottom",legend.justification ="left",panel.spacing.y =unit(0.5, "lines"),strip.placement ="outside",strip.background =element_blank() ) +guides(fill =guide_legend(nrow =1)) +scale_fill_manual(labels = new_labels_eng, values =c("#E2725B", "#556B2F", "#D2B48C", "#FFFFFF", "#808080", "#000000")) +scale_color_manual(labels = new_labels_eng, values =c("#E2725B", "#556B2F", "#D2B48C", "#FFFFFF", "#808080", "#000000"))ggsave(filename="_figs/clusters_pam_om4_q_mod_eng_25.png", seq_plot_pam_om4_q_eng, width =9.5, height =5.5, dpi=1000)
Trayectorias de hospitalización, orden de sujetos según el primer estado observado y su duración, representando a cada individuo como una línea en el gráfico (observaciones ordenadas de acuerdo a ASW)
Tiempo que demora esta sección: 0.2 minutos
Código
seq_plot2_pam_om4_q <-ggseqdplot(States_Wide.seq_quarter_t_prim_adm_cens, group= ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$clus_pam_om4,facet_ncol=2, facet_nrow=2) +theme(legend.position ="none")+# Colocar la leyenda abajolabs(x="Trimestres", y="Frecuencia relativa de estados")+theme(panel.spacing =unit(0.1, "lines"),axis.text.y =element_text(size =15), # Tamaño de las etiquetas de los grupos étnicosaxis.text.x =element_text(size =15), # Tamaño de las etiquetas del eje Xaxis.title.x =element_text(size =15), # Tamaño del título del eje Xaxis.title.y =element_text(size =15, margin =margin(r =-5)),strip.text =element_text(size =15),panel.spacing.y =unit(0.5, "lines"),strip.placement ="outside", # Para colocar las tiras fuera de los ejesstrip.background =element_blank() # Elimina el fondo para que parezca más espacioso#legend.key.size = unit(1.5, "lines"), # Aumenta el tamaño de los símbolos en la leyenda ) # Colocar la leyenda abajoseq_plot2_pam_om4_qggsave("_figs/clusterspam_om42_q_mod_25.png",seq_plot2_pam_om4_q, width =8.5, height =5.5, dpi=1000)table_data_pam_om4_q <-sprintf("%1.2f",pamRange_quarter_om$stats[3,])table_data_pam_om4_q <-as.data.frame(t(table_data_pam_om4_q))colnames(table_data_pam_om4_q)<-attr(pamRange_quarter_om$stats, "name")table_data_pam_om4_q |> knitr::kable()
PBC
HG
HGSD
ASW
ASWw
CH
R2
CHsq
R2sq
HC
0.49
0.66
0.65
0.58
0.59
948.53
0.30
1312.64
0.37
0.20
Trayectorias de hospitalización, frecuencia relativa de estados en un gráfico de barras apiladas por trimestre.
Trayectorias de hospitalización, frecuencia relativa de estados en un gráfico de barras apiladas por trimestre.
De este modo, presenta el cambio agregado en la distribución de estados a lo largo del tiempo, sin considerar las secuencias individuales.
Código
cat("Definimos las observaciones que tienen siluetas negativas")sil_neg_pam_om_clus4_q <-which(sil_pam_om_clus4_q<0)cat("A qué conglomerados pertenecen?")table(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens[sil_neg_pam_om_clus4_q, "clus_pam_om4"])ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$rn<-1:nrow(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens)
Definimos las observaciones que tienen siluetas negativasA qué conglomerados pertenecen?clus_pam_om4
6623, Un trimestre, TSM(4) 6612, Un trimestre, TUS(3)
2 0
6522, Un semestre TSM(1) 6574, Comorbilidad un trimestre(2)
151 0
Tiempo que demora esta sección: 0 minutos
1.1.2.Exploración transiciones
1.1.2.a Transiciones- RM y no RM
Tasas de transición no RM a RM y viceversa
Código
invisible("Tasas de transición no RM a RM y viceversa")trim_tasa_pam_om4_q_cens_cnt<-seqcount_t(States_Wide.seq_quarter_t_prim_adm_RM_cens, group=ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$clus_pam_om4) |> dplyr::filter(count>0) |> dplyr::mutate(trans =paste0(from,"_", to)) |> dplyr::mutate(across(c("from","to"),~gsub("\\[->\\s*|\\s*->\\s*\\]|\\[|\\]", "", .))) trim_tasa_pam_om4_q_cens_rate<-seqtrate_t(States_Wide.seq_quarter_t_prim_adm_RM_cens, group=ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$clus_pam_om4) |> dplyr::filter(rate>0) |> dplyr::mutate(trans =paste0(from,"_", to)) |> dplyr::mutate(across(c("from","to"),~gsub("\\[->\\s*|\\s*->\\s*\\]|\\[|\\]", "", .)))cat("Número de transiciones de RM a noRM y viceversa\n")cat("Transiciones de noRM a RM\n")trim_tasa_pam_om4_q_cens_cnt |>filter(from=="noRM", to=="RM") |>summarise(n=sum(count))cat("Transiciones de RM a noRM\n")trim_tasa_pam_om4_q_cens_cnt |>filter(from=="RM", to=="noRM") |>summarise(n=sum(count))cat("Total de transiciones\n")trim_tasa_pam_om4_q_cens_cnt |>summarise(n=sum(count))cat("Proporción de transiciones de noRM a RM y viceversa\n")sum(trim_tasa_pam_om4_q_cens_cnt |>filter(from=="noRM", to=="RM") |>summarise(n=sum(count)),trim_tasa_pam_om4_q_cens_cnt |>filter(from=="RM", to=="noRM") |>summarise(n=sum(count)))/trim_tasa_pam_om4_q_cens_cnt |>summarise(n=sum(count))
Tiempo que demora esta sección: 0 minutos
Código
trim_tasa_pam_om4_q_cens_rate |> dplyr::left_join(trim_tasa_pam_om4_q_cens_cnt, by=c("from"="from", "glosa_sexo"="glosa_sexo","to"="to")) |> dplyr::rename("recuento"="count") |> dplyr::filter(from %in%c("RM", "noRM")) |>ggplot(aes(x = from, y = to, fill = rate, size=log(recuento+1))) +geom_tile() +coord_flip()+scale_fill_gradient(low ="white", high ="blue") +# Ajusta la escala de colores según tus preferenciaslabs(title ="Tasas de transición, Trimestre (s/censura)",x ="Desde",y ="Hacia",fill ="Rate") +theme_minimal() +facet_wrap(~glosa_sexo)+theme(axis.text.x =element_text(angle =45, hjust =1))+geom_text(aes(label =sprintf("%1.2f", rate), size =log(recuento+1)*.5), color ="black")invisible("Hay muy pocos casos que se entrecruzan entre noRM y RM (fuera de la diagnonal)")
Porcentajes de transición no-RM y RM por cada cluster
Tiempo que demora esta sección: 0 minutos
Hay muy pocos casos que se entrecruzan entre noRM y RM (fuera de la diagnonal)
trim_tasa2_pam_om4_q_cens_rate |> dplyr::left_join(trim_tasa2_pam_om4_q_cens_cnt, by=c("from"="from", "glosa_sexo"="glosa_sexo","to"="to")) |> dplyr::rename("recuento"="count") |>#dplyr::filter(from %in% c("RM", "noRM")) |> ggplot(aes(x = from, y = to, fill = rate, size=log(recuento+1))) +geom_tile() +coord_flip()+scale_fill_gradient(low ="white", high ="blue") +# Ajusta la escala de colores según tus preferenciaslabs(title ="Tasas de transición, Trimestre (s/censura)",x ="Desde",y ="Hacia",fill ="Rate") +theme_minimal() +facet_wrap(~glosa_sexo)+theme(axis.text.x =element_text(angle =45, hjust =1))+geom_text(aes(label =sprintf("%1.2f", rate), size =log(recuento+1)*.5), color ="black")
Porcentajes de transición, transiciones posteriores, por cada cluster
Tiempo que demora esta sección: 0 minutos
1.1.2.c Tiempo promedio por cluster
Código
seq_mean_t(States_Wide.seq_quarter_t_prim_adm_cens, group=ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$clus_pam_om4) |> data.table::as.data.table(keep.rowname=T) |> dplyr::mutate(rn=gsub("\\d", "", rn)) |>ggplot(aes(x=rn, fill= factor_inclusivo, y=Mean))+geom_bar(width =1, stat ="identity") +theme_minimal() +facet_wrap(~factor_inclusivo)+labs(title =NULL,x =NULL,y =NULL) +scale_fill_manual(values =c("#70809090", "#6B8E2380", "#E27A5B","#D2B48C")) +coord_flip()+theme(#axis.text.x = element_blank(),#axis.text.y = element_blank(),panel.grid =element_blank()) +# scale_fill_brewer(palette = "Pastel1", labels=c("Sin\nautoidentificación\nni reconocimiento", "Autoidentificación\nsin reconocimiento", "Ambas")) +geom_text(aes(label =round(Mean,1)), position =position_stack(vjust =0.5), size =3.5, # Ajusta el tamaño de la fuente aquícolor ="black", # Color del textofamily ="sans", # Puedes cambiar la fuente si lo deseasbackground =element_rect(fill ="white", color =NA)) +# Fondo blancotheme(legend.title =element_blank())invisible("No me aporta mucho")
Tiempo promedio en cada estado por estatus PPOO (Trimestral c/censura)
Tiempo que demora esta sección: 0 minutos
Observamos que aquellos en el conglomerado que se encuentra un trimestre, la duración de los ingresos relacionados con trastornos de salud mental son en promedio más largos.
1.1.3. Comparación variables
1.1.3.a. Comparación covariables- PPOO
Código
ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$fac_incl_real_hist_autoperc25<-with(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens, paste0(ci_cya_conadi, ppoo_minsal_y_rsh_2010))ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens|>#dplyr::count(clus_pam_om4, factor_inclusivo_real_hist_mas_autperc)|>#2025-04-14: converted dplyr::count(clus_pam_om4, fac_incl_real_hist_autoperc25)|> dplyr::group_by(clus_pam_om4)|> dplyr::mutate(n_prop =paste0(n, " (",scales::percent(n /sum(n), accuracy=.1),")"))|> dplyr::select(-n)|> tidyr::pivot_wider(names_from = clus_pam_om4, values_from = n_prop, values_fill ="0")|> dplyr::mutate(fac_incl_real_hist_autoperc25 =factor(fac_incl_real_hist_autoperc25, levels =c("11", "10", "01", "00"), labels=c("Hay reconocimeinto/se identifica", "Hay reconocimiento/no se identifica", "No hay reconocimiento/se identifica", "No hay reconocimiento/no se identifica")))|> dplyr::arrange(fac_incl_real_hist_autoperc25)|> (\(df) {if (interactive()) {mutate(df, across(-fac_incl_real_hist_autoperc25, ~gsub("%", "", gsub("\\.", ",", .))))|> rio::export("clipboard")} knitr::kable(df, caption="Porcentajes por fila, conglomerado vs. Pertenencia/identificación + Reconocimento CONADI PPOO") })()
Porcentajes por fila, conglomerado vs. Pertenencia/identificación + Reconocimento CONADI PPOO
fac_incl_real_hist_autoperc25
6623, Un trimestre, TSM(4)
6612, Un trimestre, TUS(3)
6522, Un semestre TSM(1)
6574, Comorbilidad un trimestre(2)
Hay reconocimeinto/se identifica
336 (6.4%)
59 (7.9%)
18 (4.5%)
12 (5.2%)
Hay reconocimiento/no se identifica
107 (2.0%)
16 (2.1%)
9 (2.3%)
5 (2.2%)
No hay reconocimiento/se identifica
596 (11.4%)
92 (12.3%)
45 (11.3%)
21 (9.2%)
No hay reconocimiento/no se identifica
4209 (80.2%)
583 (77.7%)
327 (82.0%)
191 (83.4%)
Tiempo que demora esta sección: 0 minutos
Vemos las categorías de clasificación de PPOO según autopercepción (en MINSAL y en RSH) y reconocimiento CONADI.
Porcentajes por fila, conglomerado vs. Pertenencia/identificación + Reconocimento CONADI PPOO
Conglomerados
No se identifica/no pertenece
Hay reconocimeinto
6623, Un trimestre, TSM(4)
4209 (80.2%)
1039 (19.8%)
6612, Un trimestre, TUS(3)
583 (77.7%)
167 (22.3%)
6522, Un semestre TSM(1)
327 (82.0%)
72 (18.0%)
6574, Comorbilidad un trimestre(2)
191 (83.4%)
38 (16.6%)
Tiempo que demora esta sección: 0 minutos
1.1.3.b. Comparación covariables- Mortalidad
Código
# invisible("No hay nada, el tiempo promedio de censura es similar")ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens |> dplyr::mutate(death_time_rec=ifelse(death_time==20,0,1)) |> janitor::tabyl(clus_pam_om6,death_time_rec) |> dplyr::mutate(`1`=paste0(`1`," (", scales::percent(`1`/(`0`+`1`), accuracy=.1),")")) |> dplyr::left_join(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens |> dplyr::group_by(clus_pam_om6) |> dplyr::summarise(mean=sprintf("%1.1f",mean(cens_time))), by="clus_pam_om6") |> dplyr::select(-`0`)|> (\(df) {if (interactive()) {mutate(df, across(-1, ~gsub("%", "", gsub("\\.", ",", .))))|>select(2)|>t()|> rio::export("clipboard")} knitr::kable(df, "markdown", col.names=c("Conglomerado","Mortalidad observada", "Promedio"), caption="Post-hoc, conglomerado vs. Mortalidad y tiempo a censura") })()
Post-hoc, conglomerado vs. Mortalidad y tiempo a censura
Conglomerado
Mortalidad observada
Promedio
6623, Un trimestre, TSM(5)
58 (1.2%)
17.9
6612, Un trimestre, TUS(4)
19 (2.5%)
18.2
6522, Un semestre TSM(2)
8 (2.2%)
17.9
6624, TSM, 1 año después, otras causas(6)
5 (2.1%)
17.8
6574, Comorbilidad un trimestre(3)
6 (2.7%)
18.1
6268, TSM, 1 año después, TSM(1)
3 (1.7%)
18.1
Tiempo que demora esta sección: 0 minutos
Código
ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens |> dplyr::mutate(death_time_rec=ifelse(death_time==20,0,1)) |> janitor::tabyl(death_time_rec,clus_pam_om4) |> janitor::chisq.test(correct=T)#X-squared = 11.377, df = 3, p-value = 0.009854chisq_cramerv(with(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens|> dplyr::mutate(death_time_rec=ifelse(death_time==20,0,1)), table(death_time_rec , clus_pam_om4)))fisher.test(with(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens|> dplyr::mutate(death_time_rec=ifelse(death_time==20,0,1)), table(death_time_rec , clus_pam_om4), simulate.p.value=T, B=1e5))#p-value = 0.008208message("Descartando valores negativos en sil width")chisq_cramerv(with(subset(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens, !rn %in% sil_neg_pam_om_clus4_q)|> dplyr::mutate(death_time_rec=ifelse(death_time==20,0,1)), table(death_time_rec , clus_pam_om4)))# $chisq_statistic# [1] "11.31"# # $chisq_df# df # 3 # # $chisq_p_value# [1] "0.0101"# # $cramers_v# [1] "0.04"fisher.test(with(subset(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens, !rn %in% sil_neg_pam_om_clus4_q)|> dplyr::mutate(death_time_rec=ifelse(death_time==20,0,1)), table(death_time_rec , clus_pam_om4), simulate.p.value=T, B=1e5))#p-value = 0.007751invisible("no se basa en la distribución chi-cuadrado. Fisher se basa en permutaciones exactas, por lo que no se calculan df.")#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_##_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_tab_cl_mortalidad_pam_om4_q<- ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens |> dplyr::mutate(death_time_rec=ifelse(death_time==20,0,1)) |> janitor::tabyl(death_time_rec,clus_pam_om4) |>as.matrix(ncol=2)labels_pam_om4_q <-c("6623, Un trimestre, TSM(4)","6612, Un trimestre, TUS(3)","6522, Un semestre TSM(1)","6574, Comorbilidad un trimestre(2)")chisq.posthoc.test(with(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens|> dplyr::mutate(death_time_rec=ifelse(death_time==20,0,1)), table(death_time_rec , clus_pam_om4), simulate.p.value=T, B=1e5))# Realizar el análisis y crear la tabla directamentepairwise.prop.test(t(tab_cl_mortalidad_pam_om4_q[,2:5]), p.adjust.method ="holm")$p.value|>as.table() |>as.data.frame() |>rename(Grupo_1 = Var1, Grupo_2 = Var2, p_value = Freq) |>filter(!is.na(p_value)) |>mutate(Grupo_1 = labels_pam_om4_q[as.numeric(Grupo_1)],Grupo_2 = labels_pam_om4_q[as.numeric(Grupo_2)],p_value =ifelse(p_value <.001, "<.001", sprintf("%1.3f",p_value)) ) |>kable(col.names =c("Grupo 1", "Grupo 2", "Valor p ajustado"),align ="l",caption="Corrección parcial por comparaciones múltiples (Holm–Bonferroni)" )
Pearson's Chi-squared test
data: janitor::tabyl(dplyr::mutate(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens, death_time_rec = ifelse(death_time == 20, 0, 1)), death_time_rec, clus_pam_om4)
X-squared = 11.377, df = 3, p-value = 0.009854
$chisq_statistic
[1] "11.38"
$chisq_df
df
3
$chisq_p_value
[1] "0.0099"
$cramers_v
[1] "0.04"
Fisher's Exact Test for Count Data
data: with(dplyr::mutate(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens, death_time_rec = ifelse(death_time == 20, 0, 1)), table(death_time_rec, clus_pam_om4), simulate.p.value = T, B = 1e+05)
p-value = 0.008208
alternative hypothesis: two.sided
$chisq_statistic
[1] "11.31"
$chisq_df
df
3
$chisq_p_value
[1] "0.0101"
$cramers_v
[1] "0.04"
Fisher's Exact Test for Count Data
data: with(dplyr::mutate(subset(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens, !rn %in% sil_neg_pam_om_clus4_q), death_time_rec = ifelse(death_time == 20, 0, 1)), table(death_time_rec, clus_pam_om4), simulate.p.value = T, B = 1e+05)
p-value = 0.007751
alternative hypothesis: two.sided
Dimension Value 6623, Un trimestre, TSM(4) 6612, Un trimestre, TUS(3)
1 0 Residuals 3.34615360812403 -2.491149
2 0 p values 0.0066* 0.101900
3 1 Residuals -3.346153608124 2.491149
4 1 p values 0.0066* 0.101900
6522, Un semestre TSM(1) 6574, Comorbilidad un trimestre(2)
1 -1.293403 -1.429422
2 1.000000 1.000000
3 1.293403 1.429422
4 1.000000 1.000000
Corrección parcial por comparaciones múltiples (Holm–Bonferroni)
Grupo 1
Grupo 2
Valor p ajustado
6623, Un trimestre, TSM(4)
6623, Un trimestre, TSM(4)
0.047
6612, Un trimestre, TUS(3)
6623, Un trimestre, TSM(4)
0.654
6522, Un semestre TSM(1)
6623, Un trimestre, TSM(4)
0.654
6612, Un trimestre, TUS(3)
6612, Un trimestre, TUS(3)
1.000
6522, Un semestre TSM(1)
6612, Un trimestre, TUS(3)
1.000
6522, Un semestre TSM(1)
6522, Un semestre TSM(1)
1.000
Tiempo que demora esta sección: 0 minutos
Código
# Cargar las librerías necesariaslibrary(survival)library(ggplot2)# Crear la variable de supervivenciasurv_obj_4c <-Surv(time= ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$death_time,event =ifelse(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$death_time==20,0,1))# Realizar el análisis de Log-Rank (survdiff)surv_diff_4c <-survdiff(surv_obj_4c ~ clus_pam_om4,data = ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens)cat("Diferencia multiplicativa entre cluster TSM 1 trimestre y el resto")
Diferencia multiplicativa entre cluster TSM 1 trimestre y el resto
Pairwise comparisons using Log-Rank test
data: subset(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens, and clus_pam_om4 !rn %in% sil_neg_pam_om_clus4_q) and clus_pam_om4
6623, Un trimestre, TSM(4)
6612, Un trimestre, TUS(3) 0.028
6522, Un semestre TSM(1) 0.438
6574, Comorbilidad un trimestre(2) 0.348
6612, Un trimestre, TUS(3)
6612, Un trimestre, TUS(3) -
6522, Un semestre TSM(1) 1.000
6574, Comorbilidad un trimestre(2) 1.000
6522, Un semestre TSM(1)
6612, Un trimestre, TUS(3) -
6522, Un semestre TSM(1) -
6574, Comorbilidad un trimestre(2) 1.000
P value adjustment method: holm
# Crear el gráfico de Kaplan-Meier con ggplot2ggplot(km_data_4c, aes(x = time, y = surv, color = strata)) +geom_step(size =1.2) +# Curvas de supervivencia#geom_ribbon(aes(ymin = lower, ymax = upper, fill = strata), alpha = 0.2, color = NA) + # Intervalos de confianzalabs(x ="Tiempo (meses)",y ="Probabilidad de Supervivencia",color ="Grupo",fill ="Grupo" ) +theme_minimal() +ylim(c(0.9,1))+theme(legend.position ="bottom")
Tiempo que demora esta sección: 0 minutos
Código
def_enc17_21 <- readr::read_delim("_output/def_enc.csv", delim =";", escape_double =FALSE, trim_ws =TRUE)# DIAG1= Diagnóstico principal del usuario# DIAG2= Causa externa de hospitalizacióndef_enc17_21 <-def_enc17_21[,c("RUN", "DIAG1", "DIAG2")]# Suicidio (X60-X84)# Homicidio (X85-X99, Y00-Y09)# Intención indeterminada (Y10-Y34)# Armas de fuego (W32-W34)# Legal intervention (Y35.0-Y35.4, Y35.6-Y35.7)# Terrorism (U01-U03# Specific Y codes (Y87.0, Y87.1, Y87.2, Y89.9, Y86, Y89.0)invisible(paste0(" Fowler, K. A., Jack, S. P., Lyons, B. H., Betz, C. J., & Petrosky, E. (2018). Surveillance for Violent Deaths —National Violent Death Reporting System, 18 States, 2014. MMWR Surveillance Summaries, 67(2), 1. https://doi.org/10.15585/mmwr.ss6702a1"))pattern <-"^(X(6[0-9]|7[0-9]|8[0-4])\\d|X(8[5-9]|9[0-9])\\d|Y0[0-9]\\d|Y(1[0-9]|2[0-9]|3[0-4])\\d|W3[2-4]\\d|Y35[0-46-7]|Y87[012]\\d|Y89[09]\\d|Y86$|U0[1-3]\\d)"def_enc17_21 |> dplyr::filter(RUN %in%subset(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens, filter= death_time_rec==1)$run) |>#dplyr::filter(!RUN %in% runs_6025$run) |> dplyr::select(RUN, DIAG1, DIAG2) |> reshape2::melt(id.vars="RUN") |> dplyr::filter(!is.na(value)) |># janitor::tabyl(value) |> # arrange(desc(n)) |> dplyr::mutate(violent_death =grepl(pattern, value, ignore.case = F)) |> janitor::tabyl(violent_death) # violent_death n percent# FALSE 122 0.7820513# TRUE 34 0.2179487#clus_pam_om4 tstop event rate lower upper# 6623, Un trimestre, TSM(4) 104000.760 65 6.25 4.82 7.97# 6612, Un trimestre, TUS(3) 14730.962 19 12.90 7.77 20.14# 6522, Un semestre TSM(1) 7885.631 9 11.41 5.22 21.67# 6574, Comorbilidad un trimestre(2) 4487.281 6 13.37 4.91 29.10# # Fuente: https://www.suseso.cl/612/articles-18722_archivo_01.pdfdef_enc17_21 |> dplyr::filter(RUN %in%subset(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens, filter= death_time_rec==1)$run) |> dplyr::select(RUN, DIAG1, DIAG2) |> dplyr::left_join(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens[, c("run", "clus_pam_om4")], by=c("RUN"="run"), multiple="first") |> dplyr::group_by(clus_pam_om4) |> dplyr::summarise(violent_death1 =sum(grepl(pattern, DIAG1, ignore.case = F)),violent_death2 =sum(grepl(pattern, DIAG2, ignore.case = F))) |>ungroup() |> dplyr::mutate(violent_death = violent_death1 + violent_death2) |> dplyr::mutate(total = ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens |>group_by(clus_pam_om4) |>count() |>pull(n), perc= violent_death/total) |> dplyr::mutate(pq=c(104000.76, 14730.962, 7885.631, 4487.281)) |> dplyr::mutate(ir= (violent_death/pq)*1000) |> knitr::kable("markdown", caption="Tasas de mortalidad por conglomerado y causa de muerte, muertes violentas") def_enc17_21 |> dplyr::filter(RUN %in%subset(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens, filter= death_time_rec==1)$run) |> dplyr::select(RUN, DIAG1, DIAG2) |> dplyr::left_join(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens[, c("run", "clus_pam_om4")], by=c("RUN"="run"), multiple="first") |> dplyr::select(RUN, clus_pam_om4, DIAG1, DIAG2) |>mutate(categoria1 =case_when(grepl("^X(6[0-9]|7[0-9]|8[0-4])", DIAG1) ~"Lesiones autoinfligidas intencionales",grepl("^V[0-9]{2}", DIAG1) ~"Accidentes de tránsito",grepl("^(W[0-9]{2}|X[0-5][0-9])", DIAG1) ~"Otras causas externas accidentales",grepl("^(X(8[5-9]|9[0-9])|Y0[0-9])", DIAG1) ~"Agresiones",grepl("^Y8[5-9]", DIAG1) ~"Secuelas de causas externas",grepl("^Y([4-7][0-9]|8[0-4])", DIAG1) ~"Complicaciones atención médica",TRUE~"Otras causas no especificadas" ) ) |>mutate(categoria2 =case_when(grepl("^X(6[0-9]|7[0-9]|8[0-4])", DIAG2) ~"Lesiones autoinfligidas intencionales",grepl("^V[0-9]{2}", DIAG2) ~"Accidentes de tránsito",grepl("^(W[0-9]{2}|X[0-5][0-9])", DIAG2) ~"Otras causas externas accidentales",grepl("^(X(8[5-9]|9[0-9])|Y0[0-9])", DIAG2) ~"Agresiones",grepl("^Y8[5-9]", DIAG2) ~"Secuelas de causas externas",grepl("^Y([4-7][0-9]|8[0-4])", DIAG2) ~"Complicaciones atención médica",TRUE~"Otras causas no especificadas" ) ) |>group_by(clus_pam_om4) |>summarise(suicidios =sum(categoria1 =="Lesiones autoinfligidas intencionales"| categoria2 =="Lesiones autoinfligidas intencionales", na.rm =TRUE),accidentes_transito =sum(categoria1 =="Accidentes de tránsito"| categoria2 =="Accidentes de tránsito", na.rm =TRUE),otras_causas_accidentales =sum(categoria1 =="Otras causas externas accidentales"| categoria2 =="Otras causas externas accidentales", na.rm =TRUE),agresiones =sum(categoria1 =="Agresiones"| categoria2 =="Agresiones", na.rm =TRUE),secuelas =sum(categoria1 =="Secuelas de causas externas"| categoria2 =="Secuelas de causas externas", na.rm =TRUE),complicaciones =sum(categoria1 =="Complicaciones atención médica"| categoria2 =="Complicaciones atención médica", na.rm =TRUE),otras_no_especificadas =sum(categoria1 =="Otras causas no especificadas"| categoria2 =="Otras causas no especificadas", na.rm =TRUE),todas_menos_otras =sum( categoria1 %in%c("Lesiones autoinfligidas intencionales", "Accidentes de tránsito", "Agresiones") | categoria2 %in%c("Lesiones autoinfligidas intencionales", "Accidentes de tránsito", "Agresiones"),na.rm =TRUE ),.groups ="drop" ) |> dplyr::mutate(total = ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens |>group_by(clus_pam_om4) |>count() |>pull(n), perc= suicidios/total, perc_todas_menos_otras=todas_menos_otras/total) |> dplyr::mutate(pq=c(104000.76, 14730.962, 7885.631, 4487.281)) |> dplyr::mutate(ir= (suicidios/pq)*10000, ir_todas_menos_otras=(todas_menos_otras/pq)*10000) |> knitr::kable("markdown", caption="Tasas de mortalidad por conglomerado y causa de muerte", col.names=c("Conglomerado","Suicidio","Accidentes de tránsito","Otras causas externas accidentales","Agresiones","Secuelas de causas externas","Complicaciones atención médica","Otras causas no especificadas","Todas menos otras","Total","Proporción suicidio","Proporción todas menos otras", "pq", "IR suicidio", "IR todas menos otras"))
violent_death n percent
FALSE 122 0.7820513
TRUE 34 0.2179487
Tasas de mortalidad por conglomerado y causa de muerte, muertes violentas
clus_pam_om4
violent_death1
violent_death2
violent_death
total
perc
pq
ir
6623, Un trimestre, TSM(4)
0
24
24
5248
0.0045732
104000.760
0.2307675
6612, Un trimestre, TUS(3)
0
5
5
750
0.0066667
14730.962
0.3394211
6522, Un semestre TSM(1)
0
2
2
399
0.0050125
7885.631
0.2536259
6574, Comorbilidad un trimestre(2)
0
3
3
229
0.0131004
4487.281
0.6685563
Tasas de mortalidad por conglomerado y causa de muerte
pairwise_chisq_gof_test(table(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$clus_pam_om4,ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$codigo_region_rec_base), p.adjust.method="holm")|> knitr::kable("markdown", caption="Dependencia categórica sol. 4 conglomerados, por pares de categorías en RM")
Dependencia categórica sol. 4 conglomerados, por pares de categorías en RM
tab_cluster_region_pam_om4_q<-ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens |> dplyr::inner_join(data_long_establecimiento_2024_std[,c("ESTAB_HOMO", "codigo_region", "nivel_de_atencion", "nivel_de_complejidad")], by =c("estab_homo_base"="ESTAB_HOMO"), multiple ="first") |> janitor::tabyl(codigo_region, clus_pam_om4) |> janitor::adorn_percentages("col") |> janitor::adorn_rounding(digits =2)#colnames(tab_cluster_region_pam_om4_q)<- c("reg", "c1", "c4", "c3", "c5", "c6", "c7", "c8", "c9", "c2")cod_reg_homo_pam_om4_q<-data.frame(codigo_region =1:16,nombre_region =c("Región de Tarapacá","Región de Antofagasta","Región de Atacama","Región de Coquimbo","Región de Valparaíso","Región del Libertador General Bernardo O'Higgins","Región del Maule","Región del Biobío","Región de La Araucanía","Región de Los Lagos","Región de Aysén del General Carlos Ibáñez del Campo","Región de Magallanes y de la Antártica Chilena","Región Metropolitana de Santiago","Región de Los Ríos","Región de Arica y Parinacota","Región de Ñuble" ),stringsAsFactors =FALSE)dplyr::mutate(tab_cluster_region_pam_om4_q, promedio_fila =rowMeans(across(2:length(colnames(tab_cluster_region_pam_om4_q))))) |> dplyr::arrange(desc(promedio_fila)) |> dplyr::left_join(cod_reg_homo_pam_om4_q, by="codigo_region") |> dplyr::select(codigo_region, nombre_region, everything()) |> dplyr::select(-promedio_fila) |> dplyr::mutate_at(3:(length(colnames(tab_cluster_region_pam_om4_q))+1),~scales::percent(.))|> knitr::kable(caption="Porcentaje por región")
Porcentaje por región
codigo_region
nombre_region
6623, Un trimestre, TSM(4)
6612, Un trimestre, TUS(3)
6522, Un semestre TSM(1)
6574, Comorbilidad un trimestre(2)
13
Región Metropolitana de Santiago
45%
35%
46%
56%
5
Región de Valparaíso
8%
13%
11%
7%
10
Región de Los Lagos
6%
19%
4%
10%
8
Región del Biobío
10%
8%
9%
9%
9
Región de La Araucanía
5%
4%
5%
5%
6
Región del Libertador General Bernardo O’Higgins
4%
3%
5%
2%
7
Región del Maule
4%
3%
3%
3%
1
Región de Tarapacá
2%
2%
2%
2%
2
Región de Antofagasta
2%
1%
4%
1%
11
Región de Aysén del General Carlos Ibáñez del Campo
pairwise_chisq_gof_test(dplyr::filter(tab_clus_macrozona_pam_om4_q,macrozona!="RM")[-1], p.adjust.method="holm")|> knitr::kable("markdown", caption="Dependencia categórica sol. 4 conglomerados, por pares de categorías en Macrozona (corrección Holm-Bonferroni)")#Groups sharing a letter are not significantlt different (alpha = 0.05)
Dependencia categórica sol. 4 conglomerados, por pares de categorías en Macrozona (corrección Holm-Bonferroni)
pairwise_chisq_gof_test(tab_clus_sexo_pam_om4_q[-1], p.adjust.method="holm")|> dplyr::mutate(p.adj= dplyr::case_when(p.adj<0.001~"<0.001",T~sprintf("%1.3f",p.adj)))|> knitr::kable("markdown", caption="Dependencia categórica sol. 4 conglomerados, por pares de categorías en Sexo (corrección Holm-Bonferroni)")
Dependencia categórica sol. 4 conglomerados, por pares de categorías en Sexo (corrección Holm-Bonferroni)
Descriptivos, edad minima de ingreso por conglomerado
clus_pam_om4
sum
6623, Un trimestre, TSM(4)
20.58 ± 4.32 ; 20 [17, 24]
6612, Un trimestre, TUS(3)
22.41 ± 4.30 ; 23 [19, 26]
6522, Un semestre TSM(1)
19.88 ± 4.28 ; 19 [16, 23]
6574, Comorbilidad un trimestre(2)
21.78 ± 4.21 ; 21 [18, 25]
Tiempo que demora esta sección: 0 minutos
Código
dt_ing_calendar_quarter_t_desde_primera_adm_dedup|>filter(quarter ==0)|>inner_join(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens[, c("run", "clus_pam_om4")],by ="run")|> dplyr::mutate(clus_pam_om4=str_wrap(clus_pam_om4, width =20))|>group_by(clus_pam_om4)|> dplyr::summarise(mean_edad =mean(min_edad_anos),sd_edad =sd(min_edad_anos),n =n())|> dplyr::mutate(sem = sd_edad /sqrt(n),error_lower = mean_edad - sem,error_upper = mean_edad + sem)|>ggplot(aes(x = clus_pam_om4, y = mean_edad)) +geom_point(color ="black", size =3) +geom_errorbar(aes(ymin = error_lower, ymax = error_upper), width =0.2, color ="black") +labs(x ="Conglomerado (k = 4)", y ="Edad al primer evento 2018 promedio") +theme_minimal() +coord_flip() +theme(axis.text.y =element_text(size =17* .7, face ="bold"),axis.text.x =element_text(size =17* .7, face ="bold", lineheight =0.9),axis.title.x =element_text(size =16* .7, face ="bold"),axis.title.y =element_text(size =16* .7, face ="bold"),plot.margin =unit(c(0.5, 0.5, 0.5, 0.5), "cm"))ggsave("_figs/edad_minima_por_cluster_pam_om4_q_25.png", width=8, height=5, dpi=600)
Edad promedio primer ingreso con errores estándar por conglomerado
Tiempo que demora esta sección: 0 minutos
Código
invisible("Prueba de Levene par igualdad de varianzas")with(dt_ing_calendar_quarter_t_desde_primera_adm_dedup |> dplyr::filter(quarter ==0) |> dplyr::inner_join(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens[,c("run","clus_pam_om4")], by="run"), car::leveneTest(min_edad_anos, clus_pam_om4))#NSanova_clus_pam_om4_q <-oneway.test(min_edad_anos ~ clus_pam_om4, data = dt_ing_calendar_quarter_t_desde_primera_adm_dedup |> dplyr::filter(quarter ==0) |> dplyr::inner_join(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens[,c("run","clus_pam_om4")], by="run"),var.equal = T)with(dt_ing_calendar_quarter_t_desde_primera_adm_dedup |> dplyr::filter(quarter ==0) |> dplyr::inner_join(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens[,c("run","clus_pam_om4")], by="run"), oneway_anova_effect_size (min_edad_anos, clus_pam_om4))# Ver los resultados del ANOVAprint(anova_clus_pam_om4_q)#F = 49.148, num df = 3, denom df = 6622, p-value < 2.2e-16#$eta_squared#[1] 0.0217808message("Descartando valores negativos en sil width")with(dt_ing_calendar_quarter_t_desde_primera_adm_dedup |> dplyr::filter(quarter ==0) |> dplyr::inner_join(subset(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens, !rn %in% sil_neg_pam_om_clus4_q)[,c("run","clus_pam_om4")], by="run"), oneway.test(min_edad_anos ~ clus_pam_om4,var.equal = T) )#F = 48.252, num df = 3, denom df = 6469, p-value < 2.2e-16
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 3 0.4775 0.698
6622
$anova_summary
Df Sum Sq Mean Sq F value Pr(>F)
group 3 2746 915.4 49.15 <2e-16 ***
Residuals 6622 123343 18.6
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
$eta_squared
[1] 0.0217808
One-way analysis of means
data: min_edad_anos and clus_pam_om4
F = 49.148, num df = 3, denom df = 6622, p-value < 2.2e-16
One-way analysis of means
data: min_edad_anos and clus_pam_om4
F = 48.252, num df = 3, denom df = 6469, p-value < 2.2e-16
Pearson's Chi-squared test
data: df
X-squared = 30.37, df = 12, p-value = 0.002456
Fisher's Exact Test for Count Data with simulated p-value (based on
1e+05 replicates)
data: df
p-value = 0.00129
alternative hypothesis: two.sided
Porcentajes por columna, conglomerado vs. Previsión
Dependencia categórica sol. 4 conglomerados, por pares de categorías en Niv. Complejidad (corrección Holm-Bonferroni)
n
group1
group2
statistic
p
df
p.adj
p.adj.signif
5998
6623, Un trimestre, TSM(4)
6612, Un trimestre, TUS(3)
41.99007
0.00e+00
4
1.00e-07
****
5647
6623, Un trimestre, TSM(4)
6522, Un semestre TSM(1)
22.93875
1.30e-04
4
3.90e-04
***
5477
6623, Un trimestre, TSM(4)
6574, Comorbilidad un trimestre(2)
17.83641
1.33e-03
4
1.96e-03
**
1149
6612, Un trimestre, TUS(3)
6522, Un semestre TSM(1)
66.00124
0.00e+00
4
0.00e+00
****
979
6612, Un trimestre, TUS(3)
6574, Comorbilidad un trimestre(2)
18.51223
9.80e-04
4
1.96e-03
**
628
6522, Un semestre TSM(1)
6574, Comorbilidad un trimestre(2)
34.79469
5.00e-07
4
2.00e-06
****
Tiempo que demora esta sección: 0 minutos
1.1.4. Compilación comparación covariables
Código
# Definir los datos correctamentedata_pam_om4_q <-cbind.data.frame(Grupo=c("6035,\nUn trimestre,\nTSM(4)", "6025,\nUn trimestre,\nTUS(3)", "5939,\nUn semestre\nTSM(1)", "5989,\nComorbilidad un\ntrimestre(2)"), PPOO_bin =c(NA, NA, NA, NA), PPOO_sinautoid =c(NA, NA, NA, NA), PPOO_conautoid =c(NA, NA, NA, NA), Mortalidad =c(NA, NA, NA, NA), RM =c(NA, "-", NA, "+"), `Macrozona-Austral`=c(NA, NA, NA, NA), `Macrozona-Centro`=c(NA, NA, NA, NA), `Macrozona-Centro Sur`=c(NA, NA, NA, NA), `Macrozona-Norte`=c("+","-", NA, NA), `Macrozona-Sur`=c(NA, "+", NA, NA), Sexo_mujeres =c(NA, "-", NA, "-"), `Edad ingreso`=c("-", "+", "-", "+"), `Previsión-FFAA`=c("+", "-", NA, NA), `Previsión-FONASA A`=c(NA, NA, NA, NA), `Previsión-FONASA BC`=c(NA, NA, NA, NA), `Previsión-FONASA D`=c(NA, NA, NA, NA), `Previsión-ISAPRE`=c(NA, NA, NA, NA), `NivComp-Baja`=c(NA, NA, NA, NA), `NivComp-Media`=c(NA, "-", NA, NA), `NivComp-Alta`=c(NA, "+", "-", "+"))## Asegurar que los nombres de las columnas sean válidos y no haya espacios en blanco# Derretir el dataframe para que sea adecuado para ggplot2data_melt_pam_om4_q <- reshape2::melt(data_pam_om4_q, id.vars ='Grupo', variable.name ='Variable', value.name ='Asociación')# Reemplazar los NA por un valor vacíodata_melt_pam_om4_q$Asociación[is.na(data_melt_pam_om4_q$Asociación)] <-"\n"# Crear el gráfico con ggplotplot_data_melt_pam_om4_q<-data_melt_pam_om4_q |> dplyr::mutate(Variable =gsub("_", " ", Variable)) |>ggplot(aes(x = Variable, y = Grupo, fill = Asociación)) +geom_tile(color ="white", size =0.8) +scale_fill_manual(values =c("+"="#556B2F", "-"="#E2725B", "\n"="white")) +labs(title =NULL, x ="Variables", y ="Conglomerado") +theme_minimal() +theme(#axis.text.x = element_text(angle = 45, hjust = 1),panel.grid =element_blank())+theme(axis.text.y =element_text(size =17*.7, face ="bold"),#,margin = margin(l = 7)), # Tamaño de las etiquetas de los grupos étnicosaxis.text.x =element_text(size =17*.7, face ="bold"), # Tamaño de las etiquetas del eje Xaxis.title.x =element_text(size =16*.7, face ="bold"),#,margin = margin(t = -15)), # Tamaño del título del eje Xaxis.title.y =element_text(size =16*.7, face ="bold"), # Tamaño del título del eje Yplot.title =NULL, # Tamaño y estilo del título del gráficolegend.title =element_text(size =17*.7, face ="bold"), # Tamaño del título de la leyendalegend.spacing.y =unit(1.2, "lines"),legend.box.spacing =unit(0.5, "lines"), # Controla el espacio entre la leyenda y el gráficolegend.margin =margin(5, 5, 5, 5), legend.key.height =unit(1, "cm"), legend.text =element_text(size =15*.7, face ="bold") # Tamaño del texto de la leyenda ) +coord_flip()plot_data_melt_pam_om4_qggsave(filename="_figs/asociaciones_pam_om4_q_25.png", plot_data_melt_pam_om4_q, width=11*.8, height=6*.8, dpi=600)
Comparación covariables con agrupamiento
Tiempo que demora esta sección: 0 minutos
Comparación con base original PPOO
Seleccionamos casos de 2018 en adelante para hacer una base comparable.
Código
data_nueva<- rio::import("_data/EH_2010_2022_Pasantes_v2_encrip.csv")# 2) Opción: cargar tidytable DESPUÉS de dplyr para "sobreponer" sus verbslibrary(dplyr)library(stringr)library(janitor)
Adjuntando el paquete: ‘janitor’
The following object is masked from ‘package:rstatix’:
make_clean_names
The following objects are masked from ‘package:stats’:
chisq.test, fisher.test
Código
library(tidytable) # cargar después
Warning: tidytable was loaded after dplyr. This can lead to most dplyr functions being overwritten by tidytable functions.
Warning: tidytable was loaded after tidyr. This can lead to most tidyr functions being overwritten by tidytable functions.
Adjuntando el paquete: ‘tidytable’
The following objects are masked from ‘package:ggpubr’:
group_by, mutate
The following objects are masked from ‘package:rstatix’:
desc, drop_na, filter, group_by, mutate, select
The following object is masked from ‘package:magrittr’:
extract
The following objects are masked from ‘package:gtsummary’:
mutate, select
The following object is masked from ‘package:Tmisc’:
%like%
The following objects are masked from ‘package:dplyr’:
tab_clas_por_run|># tu data framemutate(clas_por_run =case_when(# 1. priorizamos Mapuche tiene_mapuche ~"Mapuche",# 2. luego otros pueblos originarios distintos de Mapuche tiene_otro_po ~"Otros reportados (no mapuche)",# 3. si no hay Mapuche ni otros y se reporta “ninguno” tiene_desconocido ~"Desconocido",# 4. si no hay Mapuche ni otros y se reporta “ninguno”!tiene_mapuche &!tiene_otro_po & tiene_ninguno ~"Ninguno puro",# 4. cualquier otro caso ⇒ NATRUE~NA_character_ ) ) |>tabyl(clas_por_run, show_na =TRUE) |>mutate(percent= scales::percent(percent, accuracy=.1))
# A tidytable: 4 × 3
clas_por_run n percent
<chr> <int> <chr>
1 Desconocido 90256 7.1%
2 Mapuche 20555 1.6%
3 Ninguno puro 1161822 91.0%
4 Otros reportados (no mapuche) 4170 0.3%
Código
# clas_por_run n percent# <chr> <int> <chr> # 1 Desconocido 90256 7.1% # 2 Mapuche 20555 1.6% # 3 Ninguno puro 1161822 91.0% # 4 Otros reportados (no mapuche) 4170 0.3% ppoo_inclusivo<- rio::import("filtered_df.csv.gz", skip=1)cat("Hacemos una clasificación comparable a la versión autorreportada de PPOO MINSAL pero para todos las hospitalizaciones del 2018")
Hacemos una clasificación comparable a la versión autorreportada de PPOO MINSAL pero para todos las hospitalizaciones del 2018
Código
tab_clas_por_run2 <- data_nueva|>filter(FECHA_INGRESO_FMT_DEIS >=as.Date("2018-01-01"), FECHA_INGRESO_FMT_DEIS <as.Date("2019-01-01"))|>replace_na(list(CyA_Conadi =0,CI_Conadi =0,RSH =0 ))|>mutate(glosa_pueblo_originario_rec =case_when( GLOSA_PUEBLO_ORIGINARIO =="NINGUNO"& RUN %in% ppoo_inclusivo[,2]~"DESCONOCIDO", T~as.character(GLOSA_PUEBLO_ORIGINARIO) # evita choque factor/chr ) )|>summarise(tiene_otro_po =any(str_detect(glosa_pueblo_originario_rec, regex_otro)),tiene_mapuche =any(str_detect(glosa_pueblo_originario_rec, "MAPUCHE")),tiene_ninguno =any(glosa_pueblo_originario_rec =="NINGUNO"),tiene_desconocido =any(glosa_pueblo_originario_rec =="DESCONOCIDO"), # si dejo al final, no recojo ningun "DESCONOCIDO" .by = RUN)#|>tab_clas_por_run2|># tu data framemutate(clas_por_run =case_when(# 1. priorizamos Mapuche tiene_mapuche ~"Mapuche",# 2. luego otros pueblos originarios distintos de Mapuche tiene_otro_po ~"Otros reportados (no mapuche)",# 3. si no hay Mapuche ni otros y se reporta “ninguno” tiene_desconocido ~"Desconocido",# 4. si no hay Mapuche ni otros y se reporta “ninguno”!tiene_mapuche &!tiene_otro_po & tiene_ninguno ~"Ninguno puro",# 4. cualquier otro caso ⇒ NATRUE~NA_character_ ) ) |>tabyl(clas_por_run, show_na =TRUE)|>mutate(percent= scales::percent(percent, accuracy=.1))
# A tidytable: 4 × 3
clas_por_run n percent
<chr> <int> <chr>
1 Desconocido 173068 13.6%
2 Mapuche 20555 1.6%
3 Ninguno puro 1079010 84.5%
4 Otros reportados (no mapuche) 4170 0.3%
#X-squared = 33.56, df = 3, p-value = 2.454e-07#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:cat("Traemos la base de datos con todos los clasificados como PPOO\n")
Traemos la base de datos con todos los clasificados como PPOO
#[1] "17.5%"cat("Y si veo del 2018 en adelante con todas las hospitalizaciones en ese periodo\n")
Y si veo del 2018 en adelante con todas las hospitalizaciones en ese periodo
Código
data_nueva|>filter(FECHA_INGRESO_FMT_DEIS >=as.Date("2018-01-01"), FECHA_INGRESO_FMT_DEIS <as.Date("2019-01-01"))|> (\(df) {nrow(distinct(df, RUN))->> filt_data_nueva df })() |>filter(RUN %in%pull(filtered_df_csv))|>distinct(RUN) |> (\(df) {print(paste0("Sólo PPOO que tienen CUALQUIERA (RSH/MINSAL histórico O CONADIS): ",nrow(df)))print(paste0("Sólo PPOO que tienen CUALQUIERA (RSH/MINSAL histórico O CONADIS) %: ",scales::percent(nrow(df)/filt_data_nueva, accuracy=.1))) })()
[1] "Sólo PPOO que tienen CUALQUIERA (RSH/MINSAL histórico O CONADIS): 197793"
[1] "Sólo PPOO que tienen CUALQUIERA (RSH/MINSAL histórico O CONADIS) %: 15.5%"
[1] "Sólo RSH o autoidentificación histórico: 151089"
[1] "Sólo RSH o autoidentificación histórico %: 11.8%"
Código
# [1] "Sólo RSH o autoidentificación histórico: 151089"# [1] "Sólo RSH o autoidentificación histórico %: 11.8%"invisible("19.9% of the patients were identified as belonging to an Indigenous Peoples group in any hospital record since 2005 or in the Social household register file (17.8%) or were officially recognized as such in the Registry of Indigenous Status/Communities (8.5%).") rm(data_nueva)rm(ci_cya_conadi)rm(filtered_df_csv)
